#include<stdio.h>
#include<stdlib.h>
#include <string.h>
#include<math.h> 
#include<iostream>
#include <direct.h>
#include<ctime>
using namespace std;

class Burgers
{
	//设置访问权限：公开
public:
		int N;
		double* macCormack(double *u0, double time, double visco, int N);
	/*A Burgers 1-demension partial differential equation numerical solver with:
		- 3 inputs: viscosity ,time and iteration times N
		- viscosity = 0, 0.01, 0.05, 0.1
		- time = 0.8, 1.4, 2.0
	Initial and boundary conditions are:
		- u(x,0) = -0.5*x, x belongs to [-1, 1]
	    - u(-1, t) = 0.5, u(1, t) = -0.5*/
};

	double* Burgers::macCormack(double *u0, double time, double visco, int N)
	{
		/*
		u is a numpy array with initial condition.
			- Predictive step: us;
			- Patial of test step: up
			- Intermediate step: unew
			- Storage of all the time step: U
		*/
    	double dx = 2.0 / N,dt = 0.0001;
		double u[N+1];
		int i, j, nt = int((double)time / (double)dt);
		for(i=0;i<N+1;i++){
			u[i]=*((double*)u0+i);
		}
		double us[N+1],up[N+1],unew[N+1];

	for(i=1;i<nt;i++){
			//Predictive step
		for(j=1;j<N;j++){
			us[j] = u[j] + (visco * dt / pow(dx,2)) * (u[j+1] - 2.0 * u[j] + u[j-1]) \
			- dt / (2.0 * dx) * (pow(u[j+1],2) - pow(u[j],2));
			}
			us[0] = 0.5;
			us[N] = -0.5;

		for(j=1;j<N;j++){
			up[j] =(visco / pow(dx,2)) * (us[j+1] - 2 * us[j] + us[j-1]) \
			- 1.0 / (2.0 * dx) * (pow(us[j],2) - pow(us[j-1],2));
			}

		for(j=1;j<N;j++){
			unew[j] = 0.5 * (u[j] + us[j] + dt * up[j]);
			u[j] = unew[j];
			}
		}
		double *pos;
    	pos=&u[0];
		return pos;
	}

int main()
{ 
	clock_t startTime,endTime;
	startTime = clock();//计时
	printf("The program is running...\n");
	int N = 200;

    //从Input.txt文档输入参数
	FILE *fp = NULL;
    char *file = (char*)"input.txt";
    char line[30];
    if( (0 != access(file,R_OK|F_OK)) || (NULL==(fp=fopen(file,"r"))) )
    {
        printf("open %s failed\n",file);
        return -1;
    }

	const char *d = ",";
    char *p;
	int off=0;
	double output[9][2];
	int kkk=-1;
    while( fgets(line, 256, fp) != NULL ){
	if(off){
    	const char *d = ",";
    	char *p;
		int kk=0;
    	p = strtok(line,d);
    	while(p)
    	{
			char str[10];
			string ing=p;
			int ii;
    	for( ii=0;ii<ing.length();ii++){
        	str[ii] = ing[ii];}
    	str[ii] = '\0';
		output[kkk][kk]=atof(str);
        p=strtok(NULL,d);
		kk++;
		}
	}
	kkk++;
	off=1;
	}
    if(fp!=NULL)
    {
        fclose(fp);
    }

	//Compute Burgers equation at different viscosity
	int i,j,ii,jj;
	Burgers burgers;
	double u[N+1],x[N+1];
	for(j=0;j<N+1;j++){
			u[j] = -0.5*(-1 + (double)j * 2/N);
			x[j] = -1 + (double)j * 2/N;
		}
	double U[kkk][N+1];
	for(i=0;i<kkk;i++){
		for(j=0;j<kkk;j++){
			
		double *pos,*poss;
		pos=&u[0];
		poss = burgers.macCormack(pos,output[j][0], output[i][1], N);
		for(ii=0;ii<N+1;ii++){
			U[j][ii]=*((double*)poss+ii);
		}
		}
		//Save 2-D data to files
		std::string prefix = "./Output";
		if (_access(prefix.c_str(), 0) == -1)
			_mkdir(prefix.c_str());   
		FILE* fpout;
		char F_PATH[30];
		sprintf(F_PATH, "./Output/burgers_mu_%4.2f.dat", output[i][1]);
		fpout = fopen(F_PATH, "w+");
		fprintf(fpout, "variables=\"x\" ");
		for(int jj = 0; jj < kkk; jj++){
			fprintf(fpout, "\"t=%4.2f\" ",output[jj][0]);
		}
		fprintf(fpout,"\nzone T=\"Burgers\" i=%d F=POINT\n",N+1);
			for (int i = 0; i < N+1; i++){
				fprintf(fpout, "%15.8f\t", x[i]);
			for(int jj = 0; jj < kkk; jj++){
			fprintf(fpout, "%15.8f\t",U[jj][i]);
			}
			fprintf(fpout, "\n");
		}
		fclose(fpout);
		printf("Finished--%d in %d\n",i+1,kkk);
	}

	endTime = clock();
	printf("Duration=%7.3fs\n",(double)(endTime - startTime)/1000);
    system("pause");
    return 0;
}